# This files estimates the discontinuity for variables expected to be unaffected to children's eligibilty

load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")

## commented out for replication files 

#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("reshape2")
#  install.packages("dplyr")
#  install.packages("rdrobust")

library(xtable)
library(dplyr)
library(rmeta)
library(ggplot2)
library(rdrobust)
library(reshape2)

mock_vars <- c("par_native",
               "par_female",
               "par_degree",
               "chi_female",
               "chi_oldest",
               "par_cohabit")
year <- c(2009, 2013, 2014, 2015)

coefmat <- array(NA, dim=c (6,5,4))

data_master <- left_join(data_day_year, data_n_year)
for (j in 1:6){
  ## Create large dataset 
  ## Define treatment dummy and assign weights 
  var <- mock_vars[j]
  data_close <- 
    data_master %>%
    ungroup() %>%
    mutate(voted     = round(unlist(data_master[,var])*n),
           abstained = round((1-unlist(data_master[,var]))*n),
           treated   = days > 0) %>% 
    select(year, days, treated, voted, abstained)
  
  data_voted <- 
    as.data.frame(lapply(data_close, 
                         function(x,p) rep(x,p), 
                         data_close[["voted"]])) %>%
    mutate(voted = 1)
  
  data_abstained <- 
    as.data.frame(lapply(data_close, 
                         function(x,p) rep(x,p), 
                         data_close[["abstained"]])) %>%
    mutate(voted = 0)
  
  data <- 
    rbind(data_voted, data_abstained) 
  
  for (i in 1:4){
    year_vec <- year[i]
    data_year <- 
      data %>%
      mutate(days = days - 0.5) %>% 
      filter(year == year_vec )
    
    rd <- rdrobust(data_year$voted, data_year$days)
    coefmat[j,i,c(1,3)] <- rd$coef[c(1,3)]
    coefmat[j,i,c(2,4)] <- rd$se[c(1,3)]
    
  }
  
  coefmat[j, 5, 1] <- meta.summaries(coefmat[j, 1:4, 1], coefmat[j, 1:4, 2])$summary
  coefmat[j, 5, 2] <- meta.summaries(coefmat[j, 1:4, 1], coefmat[j, 1:4, 2])$se.summary
  coefmat[j, 5, 3] <- meta.summaries(coefmat[j, 1:4, 3], coefmat[j, 1:4, 4])$summary
  coefmat[j, 5, 4] <- meta.summaries(coefmat[j, 1:4, 3], coefmat[j, 1:4, 4])$se.summary
  print(j)
}

#bind estimates and standard error in data frame
plot_data1 <- 
  data.frame(rbind(coefmat[1, , c(1, 2)],
                   coefmat[2, , c(1, 2)],
                   coefmat[3, , c(1, 2)],
                   coefmat[4, , c(1, 2)],
                   coefmat[5, , c(1, 2)],
                   coefmat[6, , c(1, 2)],
                   coefmat[1, , c(3, 4)],
                   coefmat[2, , c(3, 4)],
                   coefmat[3, , c(3, 4)],
                   coefmat[4, , c(3, 4)],
                   coefmat[5, , c(3, 4)],
                   coefmat[6, , c(3, 4)]))

colnames(plot_data1) <- c("est", "se") # assign names to variables

# create other variable
plot_data <- 
  plot_data1 %>% 
  mutate(rd     = rep(c("Conventional RD", "Robust RD"), each = 30),
         year   = rep(c(year, "Pooled"), 12),
         alpha  = rep(c(rep(0.5, 4), 1), 12),
         place  = rep(1:5, 12),
         var    = rep(rep(c("Parent is \nnative",
                            "Parent is \nfemale",
                            "Parent has \nbachelor's",
                            "Child is \nfemale",
                            "Child is \noldest",
                            "Parent and child \ncohabit"), each = 5), 2))

plot_data$var <- factor(plot_data$var, 
                        levels = c("Parent is \nnative",
                                   "Parent is \nfemale",
                                   "Parent has \nbachelor's",
                                   "Child is \nfemale",
                                   "Child is \noldest",
                                   "Parent and child \ncohabit"))


plot <- 
  ggplot(data = plot_data,
         aes(x = 100*est, 
             xmin = 100 * (est + qnorm(0.025) * se),
             xmax = 100 * (est + qnorm(0.975) * se),
             y    = place)) + 
  geom_point(alpha = plot_data$alpha) + 
  geom_errorbarh(alpha = plot_data$alpha, height = 0) +
  facet_grid(var ~ rd) + 
  theme_classic() + 
  scale_y_continuous("", labels = c(year, "Pooled")) + 
  scale_x_continuous("Effect in percentage points") + 
  theme(panel.spacing = unit(1, "lines")) + 
  geom_vline(xintercept = 0, alpha = 0.3, linetype = 2)